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The  paper  presents  a  simple  approximation  technique  for  statistical  tolerance 
analysis,  namely,  the  allocation  of  component  tolerances  based  on  a  known 
assembly  tolerance.  The  technique  utilizes  a  discretized,  multivariate  kernel  den¬ 
sity  estimate  and  a  simple  transformation  to  approximate  the  probability  distri¬ 
bution  of  the  overall  assembly  characteristic.  The  data-driven  approach  is  suitable 
for  real-world  settings  in  which  components  are  randomly  selected  from  their 
respective  manufacturing  processes  to  form  mechanical  assemblies. 
Demonstrated  is  the  numerical  approach  in  two  dimensions  for  two  distinct 
cases:  first,  when  component  characteristics  are  non-normal,  independent 
random  variables,  and  second,  when  they  are  highly  correlated,  normal  random 
variables.  The  results  are  promising  in  initial  test  problems. 


1.  Introduction 

Many  mechanical  assemblies  consist  of  two  or  more  components  or  subassem¬ 
blies.  Owing  to  variation  in  machine  performance,  operator  performance,  or  proper¬ 
ties  of  the  work  piece,  it  is  impossible  to  avoid  variation  in  component  characteristics 
completely.  Naturally,  variation  in  the  components  leads  to  variation  in  the  overall 
assembly;  thus,  the  component  characteristics  (and  the  assembly  characteristic)  may 
be  considered  as  continuous  random  variables.  When  the  assembly  characteristic  is  a 
linear  function  of  the  component  characteristics,  normality  and  independence 
assumptions  allow  for  a  simple  calculation  of  the  component  tolerances  to  yield  a 
desirable  acceptance  rate  for  the  assemblies.  However,  in  many  situations,  neither 
the  component  nor  assembly  characteristic  probability  distributions  are  known. 
Furthermore,  underlying  dependencies  may  exist  between  component  characteristics 
that  distort  the  results  of  the  aforementioned  analysis. 

This  paper  illustrates  an  approximation  technique  for  statistically  allocating 
component  tolerances,  based  on  a  prespecified  assembly  tolerance,  when  the  prob¬ 
ability  distributions  of  component  and  assembly  characteristics  are  not  known  or 
when  the  component  characteristics  are  normally  distributed  and  correlated.  It  is 
further  assumed  that  the  assembly  characteristic  is  a  linear  function  of  the  com¬ 
ponent  characteristics  or  a  form  that  may  be  linearized.  The  next  section  reviews 
some  of  the  past  and  current  literature  in  the  area  of  statistical  tolerancing  analysis; 
section  3  presents  the  approach  to  the  problem,  while  section  4  gives  a  review  of 
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kernel  density  estimation;  section  5  presents  numerical  examples  demonstrating  the 
technique  when  the  assembly  is  comprised  of  two  components,  and  section  6  gives 
concluding  remarks  and  directions  for  future  research. 

2.  Review  of  past  work  and  problem  statement 

The  focus  of  the  work  is  on  statistical  tolerancing  analysis.  Some  other  broadly 
defined  tolerancing  techniques  include  manual  tolerancing  based  on  the  dimensional 
tolerancing  chain,  computer-aided  tolerancing  techniques  (Bjorke  1989),  and  toler- 
ancing-based  geometrical  modelling.  As  noted  by  Zhang  and  Huq  (1992),  software 
implementations  of  these  techniques  have  been  developed.  In  the  usual  way,  it  is 
assumed  that  the  overall  assembly  characteristic,  X ,  is  given  as  a  function  of  the 
component  characteristics,  i  =  1, 2, . . . ,  k  of  the  form 

X  =  ^(Xl,X2,...,Xk),  (1) 

where  X  is  assumed  to  have  density  /,  while  X ,•  has  density  /•  for  i  =  1, 2, . . . ,  k.  In 
most  analyses,  it  is  assumed  that  the  component  characteristics  are  independent 
random  variables  with  known  densities.  However,  characterization  of  the  density, 
/,  is  difficult  in  the  absence  of  these  assumptions.  For  this  reason,  techniques  for 
finding  the  moments  of  the  distribution  have  been  proposed. 

Many  authors,  such  as  Burr  (1958),  Mansoor  (1963),  Bender  (1968)  and  Evans 
(1974),  have  considered  the  problem  when  if;  is  linear.  Evans  (1974),  Nigam  and 
Turner  (1995)  and  Barker  (1989)  discussed  the  linear  propagation  of  errors  or 
stack  tolerancing  in  which  the  first  two  moments  of  the  distribution  of  X  are  deter¬ 
mined  by  partial  differentiation  of  ^  with  respect  to  each  component  characteristic. 
Mansoor  (1963)  considered  the  tolerance  problem  in  engineering  design  and 
assumed  that  the  overall  assembly  tolerance  is  prespecified  by  customer  or  engin¬ 
eering  design  requirements.  His  objective  was  to  allocate  the  component  tolerances 
so  that  the  assembly  tolerance  objective  was  met,  yielding  an  acceptable  rejection 
rate.  The  technique  resulted  in  a  probabilistic  allocation  of  component  tolerances, 
but  was  heavily  contingent  upon  the  assumption  of  normality  and  independence  of 
component  characteristics.  Parkinson  (1982)  extended  Mansoor’ s  work  by  estimat¬ 
ing  the  probability  that  the  component  characteristics  would  fail  to  meet  the  dimen¬ 
sional  requirements  of  the  assembly  characteristic.  Parkinson  (1984)  optimized 
component  tolerances  so  that,  for  a  given  manufacturing  cost,  the  probability  of 
components  failing  to  meet  the  assembly  specification  limits  was  below  some  accep¬ 
table  level.  The  technique  requires  the  component  characteristics  to  follow  either  a 
normal,  truncated  normal  or  mixed  normal  distribution. 

When  the  function  ?/>  is  non-linear,  it  is  necessary  to  employ  the  non-linear 
propagation  of  errors  or  extended  Taylor  series  approximation.  Evans  (1974) 
showed  that  if  the  derivatives  of  tp  were  tractable,  then  the  first  four  moments  of 
the  distribution  of  X  may  be  found  analytically  without  too  much  difficulty. 
However,  when  ^  is  unknown  analytically,  one  must  resort  to  numerical  integration 
to  evaluate  the  moments  of  the  distribution.  Evans  (1974)  and  Nigam  and  Turner 
(1995)  gave  the  appropriate  quadrature  expression  for  the  mean  of  X;  however,  the 
technique  assumes  independence  of  the  component  characteristics.  In  using  the 
quadrature  method,  one  advantage  is  that  the  distribution  of  component  character¬ 
istics  need  not  be  known,  only  their  lower  moments.  D'Errico  and  Zaiano  (1988) 
described  a  modification  to  Taguchi's  method  for  determining  the  moments  of  the 
distribution  of  X.  Specifically,  component  distributions  were  assumed  to  be  indepen- 
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dent  and  normal,  but  were  estimated  by  a  discrete  distribution  at  three  points.  The 
method  uses  a  modified  (Gaussian)  quadrature  approach  that  outperforms  Taguchi's 
method,  at  least  for  the  linear  case,  by  correctly  specifying  the  fifth  moment  of  the 
distribution.  More  recently,  the  concept  of  distribution  function  (DF)  zones  has  been 
introduced  into  the  statistical  tolerancing  literature  (cf.  O'Connor  and  Srinivasan 
1997).  In  this  approach,  a  population  of  components  is  deemed  acceptable  if  its 
distribution  function  is  bounded  above  and  below  by  some  specific  distribution 
functions.  This  part-level  technique  has  been  extended  by  Zhang  et  al.  (1999)  to  a 
statistical  tolerance  synthesis  model  that  additionally  minimizes  the  manufacturing 
cost.  The  DF  zone  approach  requires  that  component  characteristics  are  statistically 
independent. 

The  present  paper  contributes  to  the  literature  on  statistical  tolerancing  analysis 
by  considering  two  separate  cases  in  the  allocation  of  component  tolerances  when 
the  assembly  tolerance  is  fixed.  First,  it  considers  the  case  in  which  a  small  number  of 
component  characteristics  are  independent  and  follow  non-normal  distributions.  In 
the  second  case,  it  considers  normally  distributed  component  characteristics  that  are 
highly  correlated.  We  present  an  approximation  technique  to  allocate  the  tolerance 
of  the  assembly  characteristic  among  k  components  so  that  an  overall  desired  assem¬ 
bly  tolerance  may  be  achieved  when  ^  is  a  linear  function  of  the  k  component 
characteristics.  The  data-driven  technique  uses  a  non-parametric  estimate  of  the 
joint  density  of  component  characteristics  to  approximate  the  distribution  function 
of  X  rather  than  assuming  a  parametric  distribution. 


3.  Approximation  technique 

The  following  assumptions  are  made. 

•  The  assembly  characteristic,  X,  and  component  characteristics,  Xh 
i  =  1,2 ,k  are  of  the  Nominal-the-Better  type  (N-type). 

•  Xj  and  Xj,  i  /  /  are,  in  general,  correlated  random  variables  for  i, 
j  =  1,2  ,...,k. 

•  Components  are  randomly  assembled. 

•  Component  characteristics  follow  some  unknown  probability  distribution. 

•  The  process  that  generates  Xt  is  adjusted  so  that  the  mean  of  Xh  p,h  equals  the 
nominal  size  of  Xh  i  =  1,2 ,k. 

•  The  proportion  of  Xt  falling  within  the  specification  limits  is  1  —  /;,  for 

/=  1,2, ... 

•  The  proportion  of  X  falling  within  the  specification  limits  is  1  —  p. 

•  The  assembly  tolerance  T  is  specified  by  the  customer  or  by  engineering 
requirements. 

•  The  specification  limits  need  not  be  equidistant  from  the  nominal  dimension  of 
the  characteristic. 

The  natural  process  capability  of  the  process  that  generates  the  components  and  the 
assembly  is  a  proportion  of  the  range  of  all  possible  values  that  the  component  or 
assembly  characteristics  may  assume.  Let  /,•  denote  the  natural  process  capability  of 
the  process  that  generates  component  characteristic  Xx  and  let  t  denote  the  natural 
process  capability  for  the  process  that  generates  the  assembly  characteristic  X.  Now 
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for  any  component  or  assembly  characteristic  distribution,  the  natural  process  cap¬ 
ability  may  be  given  as  a  function  of  the  process  standard  deviation  by 

h  =  Ci<Jh  i  =  1,2, ...  ,*  (2) 

and 

t  =  ca,  (3) 

where  c,  and  c  are  constants,  07  is  the  standard  deviation  for  the  process  that  gen¬ 
erates  component  i,  and  a  is  the  standard  deviation  for  the  process  that  generates  the 
assembly  characteristic.  The  values  of  ct  and  c  dictate  the  proportion  of  the  total 
distribution  range  that  is  covered  by  the  natural  process  capability  and  depend  upon 
the  distribution  shape.  For  instance,  if  the  component  characteristics  are  normally 
distributed,  then  the  assembly  characteristic  is  also  normally  distributed  and  c  =  6 
yields  a  range  of  values  for  the  assembly  characteristic  that  covers  99.73%  of  the 
entire  distribution  range  (0.135%  on  each  tail).  If  c  =  8,  then  99.994%  of  the  distri¬ 
bution  is  covered.  Let  7)  denote  the  tolerance  of  component  i  and  let  T  denote  the 
assembly  tolerance.  The  process  capability  ratio,  denoted  by  Cp,  for  either  a  com¬ 
ponent  or  assembly  characteristic  is  the  ratio  of  the  width  of  the  tolerance  interval  to 
the  natural  process  capability  (e.g.  Cp  =  T /t).  It  is  assumed  here  that  the  manufac¬ 
turer  specifies  a  target  process  capability  ratio,  C*,  for  the  component  and  the 
assembly  characteristics.  Using  the  substitution  g  =  cC*  and  gt  =  c;C*, 
i  =  1, 2, . . . ,  k  it  is  seen  that 

T  =  U  -  L  =  ga 

Tj  =  Uj  —  Lj  =  gi<jj  1=1,2,...,*,  (4) 

where  U  and  L  are  the  upper  and  lower  specification  limits  of  the  assembly,  respect¬ 
ively,  and  Uj  and  Lj  are  the  upper  and  lower  specification  limits  for  component  i, 
respectively.  It  can  be  seen  by  (4)  that  the  component  tolerances  are  not  equal  in 
general,  but  rather,  each  is  uniquely  determined  by  g  or  g,-,  which  directly  depends  on 
the  distribution  of  the  component  or  assembly  characteristic.  The  following  subsec¬ 
tion  relates  the  prespecified  assembly  tolerance  to  the  k  component  tolerances  while 
the  means  for  computing  the  constant  g  and  g,-,  i  =  1,2, ...  ,*  will  be  presented  in 
section  3.2. 


3.1.  Relating  assembly  and  component  tolerances 

Let  p  and  a2  denote  the  mean  and  variance  of  the  assembly  characteristic  X , 
respectively.  For  any  set  of  component  characteristics,  whether  correlated  or  uncor¬ 
related,  it  is  well  known  that 

k 

M  =  X  Vi-  (5) 

1=1 

The  variance  of  the  assembly  characteristic  is  given  by  (Freund  and  Walpole  1987) 

k  k—\k 

cr  =  ^2  a2 a2  +  2^2  (6) 

i=l  i=l  7=i+l 

where  a,-  is  a  constant,  aj  is  the  variance  of  the  zth  component  characteristic  and  af  ,• 
is  the  covariance  of  Xt  and  Xj,i^  j.  In  case  the  component  characteristics  are 
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independent,  at  =  1  and  <t,-  ;  =  0  for  all  i,  j  =  1,2,...,  k.  The  fifth  assumption  of  the 
problem  asserts  that  the  process  generating  component  i  is  adjusted  so  that  the  mean 
p,j  is  equal  to  the  nominal  dimension,  say  dj.  However,  it  is  likely  that  during  process 
setting,  the  operator  may  not  be  able  to  adjust  the  process  mean  to  the  nominal 
setting  cl,.  To  account  for  this  phenomenon,  one  may  adopt  the  approach  of 
Mansoor  (1963)  who  assumed  a  ‘worst-case’  tolerance  wherein  the  zth  mean  can 
shift  no  more  than  a  fixed  distance  from  the  nominal  value.  This  mathematically 
imposes  the  conditions 

\Hi-di\<sh  i=  1,2,...,*,  (7) 

where  s ,•  denotes  the  maximum  allowable  deviation  from  the  nominal  dimension  of 
the  zth  component.  It  can  easily  be  shown  (Mansoor  1963)  that  the  assembly  toler¬ 
ance  specified  by  design  engineers  to  reflect  these  modifications  would  be  given  by 

T'=T  -2X>  (8) 

;=1 


Equation  (8),  in  effect,  would  reduce  the  amount  of  tolerance  that  may  be  allocated 
among  the  k  component  characteristics.  Now,  substituting  (6)  into  (4),  we  obtain 


T  =  g 


(9) 


Equation  (9)  provides  one  equation  relating  the  tolerance  of  the  assembly  character¬ 
istic  to  the  component  characteristics,  Th  i=  1,2 ,...,k.  Thus,  we  require  k  —  1 
additional  equations  to  solve  for  the  component  tolerances,  and  these  are  obtained 
by  the  assumption  that  the  target  process  capability  ratio,  C*,  is  equal  for  each  of  the 
processes  generating  the  k  component  characteristics.  By  (2)  and  (3)  and  the  sub¬ 
stitutions  g  =  cC*  and  g ,■  =  c,C*,  i  =  1,2, ...  ,k,  this  assumption  implies 


T\  T2  Tk 

S  lCT,  g2cr2  gkak 


(10) 


Now,  (9)  combined  with  any  k  —  1  equations  in  (10)  constitutes  a  set  of  k  equations 
that  may  be  solved  simultaneously  to  allocate  probabilistically  the  assembly  toler¬ 
ance  T  among  the  component  tolerances,  Th  i  =  1,2,...,  k.  However,  this  allocation 
must  take  into  account  the  desired  proportion  of  rejects  for  the  components  and  the 
final  assembly.  The  following  subsection  describes  the  procedure  for  finding  g  and  gh 
i=  1,2,  ...,k,  for  a  desired  proportion  of  rejects  and  for  allocating  the  assembly 
tolerance  among  the  component  tolerances. 


3.2.  Tolerance  allocation  procedure 

Let  the  cumulative  distribution  functions  (CDFs)  of  Xh  i  =  1,  2, . . . ,  k  and  X  be 
Fx,  i  =  1,2 ,k  and  Fx,  respectively.  Define  pf1  as  the  proportion  of  undersized 
components  of  type  i  and  define  py  as  the  proportion  of  oversized  components  of 
the  same  type.  The  total  proportion  of  rejects  for  the  process  that  generates  the  zth 
component  is  given  by  p ,•  =  py  +  p^\  Similarly,  denote  the  proportion  of  oversized 
assemblies  by  pa ,  and  the  proportion  of  undersized  assemblies  by  pu  such  that 
P  =  Pu  +  /V  where  p  denotes  the  total  proportion  of  rejects  for  the  assembly.  The 
upper  and  lower  specification  limits  need  not  be  equidistant  from  the  mean  and  the 
distribution  may  not  be  symmetric;  hence,  the  oversized  and  undersized  proportions 
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Figure  1.  Depiction  of  the  variables  for  an  arbitrary  distribution. 


need  not  be  equal.  Figure  1  depicts  the  variable  definitions  for  an  arbitrary  assembly 
characteristic  distribution. 

Equations  (9)  and  (10)  indicate  that  gh  i  =  1,2, . . . ,  k  and  g  are  needed  to  allocate 
the  assembly  tolerance  T  among  the  k  component  tolerances.  The  value  of  gh 
i  =  1, 2, . . . ,  k  is  obtained  from  the  cumulative  distribution  function  of  component 
characteristic  i,  namely  Fx,  i=  1,2,  ...,k,  while  g  is  obtained  from  the  assembly 
characteristic  CDF,  Fx,  as  outlined  in  the  following  procedure  description. 

Step  0.  Choose  acceptable  values  py,  py,  pu,  p„  and  C*. 

Step  1.  Compute  c  and  c,  for  the  assembly  and  each  component  characteristic  by 


where 


-  u 

fV)  _ 

A  lu  ll 

and  ct  = - — , 

(U) 

a 

sup{x 

:  Fx{x)  <  pu } 

(12) 

sup{.\'(. 

:  l'\  (.v,  l  <  /)„'  } 

(13) 

infj.v  : 

1  -  Fx(x)  <  p„ } 

(14) 

inf  {-T 

:  1  -fx,{Xt)<pf}. 

(15) 

Step  2.  Calculate  g  and  g ,•  for  the  assembly  and  each  component  characteristic  by 

g  =  cC*  and  g,  =  c,C*.  (16) 

Step  3.  Obtain  Th  i  =  1,2 ,...,k  by  simultaneously  solving  (9)  and  (10). 


The  procedure  here  assumes  the  cumulative  distribution  functions  of  components 
and  the  assembly  are  known  a  priori ;  however,  this  is  seldom  the  case  in  reality.  One 
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alternative  to  the  normality  assumption  is  to  utilize  data-driven,  non-parametric 
distribution  approximations,  and  one  well-known  technique  is  reviewed  in  section  4. 


4.  Distribution  approximations 

In  the  absence  of  normality  (and  independence)  assumptions  of  the  component 
characteristics,  it  is  difficult  to  find  an  analytical  solution  for  the  cumulative  distri¬ 
bution  function  of  the  assembly  characteristic,  X.  Thus,  numerical  approximations 
will  be  employed  to  carry  out  this  task.  In  this  section,  a  well-known,  non-parametric 
density  estimation  technique  and  its  multivariate  extension  are  reviewed  to  find  the 
distribution  of  X. 


4. 1 .  Univariate  kernel  density  estimator 

A  brief  review  of  the  kernel  density  estimate  (KDE)  is  presented  (cf.  Silverman 
1986).  Suppose  a  continuous  random  variable  Y  has  probability  density  function 
f(y)  and  let  Yh  i  =  1,2,...,  n,  denote  n  independent,  identically  distributed  observa¬ 
tions  of  Y.  The  kernel  density  estimate  of  /  at  a  point  y  is  given  by 

Z—  1  x  7 


where  K  is  the  kernel  function,  n  is  the  number  of  real  observations  and  h  is  the 
smoothing  parameter  or  bandwidth.  The  density  estimate  at  a  fixed  point  (y)  is  the 
sum  of  n  kernels  of  mass,  where  the  magnitude  of  each  kernel  is  dictated  by  the 
distance  between  the  fixed  point  and  the  observed  data  point  (i.e.  observations  that 
are  within  the  /?-radius  of  the  fixed  point  contribute  more  to  the  density  estimate  at 
that  point  while  those  outside  the  radius  contribute  little).  The  kernel  function  is 
usually  (but  not  necessarily)  a  probability  density  function  itself,  being  positive 
valued  on  the  real  line  and  integrating  to  unity.  Its  shape  dictates  the  shape  of 
each  kernel  while  the  smoothing  parameter  determines  its  width.  A  large  smoothing 
parameter  yields  a  smoother  density  estimate  while  a  smaller  value  reveals  fine  details 
in  the  density  estimate.  One  common  kernel  function  for  univariate  data  is  the 
Gaussian  kernel  given  by 


Kg{u) 


1  ^ 

■  exp  (— 0.5w  ). 

v27T 


(18) 


As  is  the  case  with  most  statistical  estimators,  there  exists  a  trade-off  between  the 
variance  and  bias  of  a  KDE.  This  trade-off  is  dictated  by  adjustment  of  the  smooth¬ 
ing  parameter  /?,  and  when  the  Gaussian  kernel  function  is  used,  the  optimal  smooth¬ 
ing  parameter  may  be  computed  directly  by  (Silverman  1986) 

/?opt  =  0.9  ■  min  {a,  interquartile  range/1.34}  •  (19) 

where  cr  may  be  approximated  with  the  sample  standard  deviation,  s.  Silverman 
(1986)  stated  that  this  value  works  well  for  a  wide  range  of  univariate  distributions, 
including  those  with  a  bimodal  structure.  Furthermore,  the  author  demonstrated 
through  simulation  that  for  n  =  100,  skewness  or  bimodality  will  be  revealed  using 
the  smoothing  parameter  of  (19).  However,  one  drawback  of  the  kernel  estimator  is 
its  tendency  to  show  spurious  noise  in  long-tailed  distributions.  This  is  due  to  the  fact 
that  a  fixed  smoothing  parameter  is  used  over  the  range  of  the  data.  By  employing 
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adaptive  kernel  methods  which  vary  the  smoothing  parameter  over  the  range  of  the 
data,  this  problem  can  be  alleviated. 

We  apply  the  univariate  KDE  of  (17)  in  the  following  manner.  Let  P  be  a  finite 
set  of  fixed  points  at  which  to  estimate  the  density.  The  probability  mass  associated 
with  the  point  x0  G  P  is  given  by 

P{*  o)  =  S~lf(x0),  (20) 

where  6  =  Sxep/M-  The  cumulative  distribution  function  approximator  is  given  by 

F(*o)  =  (21) 

Assuming  independence  of  the  component  characteristics,  it  is  possible  to  obtain  the 
distribution  of  X  by  means  of  the  convolution  of  the  component  characteristic 
density  functions,  or  by  the  central  limit  theorem  if  there  are  a  large  number  of 
components.  However,  since  the  component  characteristics  are  not  necessarily  inde¬ 
pendent,  the  distribution  must  be  obtained  by  means  of  the  joint  density.  Next,  the 
multivariate  kernel  estimator  will  be  reviewed  for  approximating  the  joint  density 
function  of  the  component  characteristics. 


4.2.  Multivariate  kernel  density  estimator 

The  distribution  of  the  assembly  characteristic  is  not  easily  obtained,  particularly 
when  the  assumption  of  independence  is  removed.  In  this  section,  the  multivariate 
form  of  the  kernel  density  estimator  is  reviewed  (Silverman  1986,  Scott  1992). 

Let  X1 ,  X2, , . . ,  Xk  be  k  continuous  random  variables  having  joint  density  func¬ 
tion/^,  x2>  . . .  Further  suppose  that  the  random  variables  are  a  function  of  X 
such  that 

X  =  ^(XuX2,...,Xk).  (22) 

The  density  function  of  X  may  be  obtained  from  the  joint  density  of  the  components 
by  the  moment-generating  function  (MGF)  technique  or  the  transformation  of  vari¬ 
ables  technique.  The  transformation  of  variables  technique  is  employed  since  the 
MGF  technique  requires  independence  of  the  random  variables.  The  multivariate 
form  of  the  kernel  estimator  provides  the  means  by  which  the  joint  density  of  com¬ 
ponent  characteristics,  f(xl,  x2,  ...,  xk)  may  be  approximated  and  is  given  by 

<*> 

where  x  is  a  point  in  /c-dimensional  space.  Thus,  x;-  =  (xy,  x2j,  . . . ,  Ay),./  =  1,2 
constitutes  a  set  of  n  observations  from  the  manufacturing  process  producing  com¬ 
ponent  Xh  i  =  1,2 ,k.  It  is  seen  that  the  kernel  estimate  for  multivariate  data  is 
analogous  to  the  univariate  case  in  that  we  consider  the  ‘distance’  from  our  point  of 
consideration,  x  and  the  observations,  x;-,  j  =  1,2 ,...,«.  Furthermore,  the  kernel 
function  is  positively  valued  and  integrates  to  unity  on  R/l. 

If  there  exist  extreme  differences  in  variance  between  the  component  character¬ 
istic  observations,  (23)  may  not  yield  a  good  estimate  for  the  overall  joint  density. 
Silverman  (1986)  suggested  a  method  for  prescaling  the  data  so  that  the  complication 
may  be  alleviated.  More  specifically,  the  author  suggested  a  transformation  of  the 
form  (Fukunaga  1972) 
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ici-i/2  n 

/w  =  -  x/)Ts_1(x  -  x/)h  (24) 

1 

where  S  represents  the  covariance  matrix  for  the  observed  data,  |S|  is  its  determinant 
and  is  its  inverse.  Using  a  substitution  of  variable,  let 

ui  =  /U2(x  —  x,)t5_1(x  —  x,-),  i  =  1, 2, . . . ,  h.  (25) 


Each  kernel  mass  is  obtained  by  evaluating  K  at  the  point  ur  A  common  kernel 
function  in  the  multivariate  case  is  the  multivariate  Epanechnikov  given  by 
Silverman  (1986): 


Ke(x)  =  { 


—  (k  +  2)(1  —  xTx)  xTx  <  1 

2  vk 

0  elsewhere, 


(26) 


where  k  is  the  number  of  variables  and  vk  is  the  volume  of  the  ^-dimensional  unit 
sphere.  For  example,  in  case  k  =  2,  the  value  of  the  kernel  function  is 


KM ) 


Ui)  ut  <  1 

*  =  1,2, ...  ,w 

otherwise. 


(27) 


Silverman  (1986)  and  Scott  (1992)  have  stated  that  selection  of  the  kernel  function  is 
not  a  crucial  issue  in  constructing  the  density  estimate,  thus,  the  Epanechnikov 
kernel  is  used  throughout  this  work  due  to  its  computational  simplicity. 

Selection  of  the  smoothing  parameter  for  the  multivariate  kernel  is  analogous  to 
that  of  the  univariate  case.  The  objective  again  is  to  optimize  the  trade-off  between 
variance  and  bias  in  the  estimator.  The  optimum  smoothing  parameter  is  given  by 

V  =  A{K)iMk+ 4>"1>  (28) 

where  A(K)  is  a  constant  that  depends  on  the  kernel  function.  For  example, 
A(K)  =  2.40  when  using  the  Epanechnikov  kernel  for  the  case  k  =  2  (cf. 
Silverman  1986).  The  link  between  the  joint  density  of  the  component  characteristics 
and  the  density  of  the  assembly  characteristic,  X,  is  made  through  the  transformation 
of  variables  technique.  Implementation  of  the  technique  requires  the  continuous 
characteristic  distributions  to  be  approximated  by  discrete  distributions;  therefore, 
the  discrete  form  of  the  transformation  should  be  applied. 


5.  Validation  and  numerical  examples 

Distributions  of  the  component  and  assembly  characteristics  will  be  estimated  on 
a  finite  support;  thus,  the  transformation  of  variables  technique  for  discrete  variables 
is  reviewed  for  the  case  k  =  2.  Let  X{  and  X2  be  two  discrete  random  variables  with 
joint  probability  mass  function  (pmf),  <f>.  The  transformation  of  variables  technique 
transforms  the  joint  pmf  of  Xx  and  X2  to  the  joint  pmf  of  X  and  Xx.  The  procedure 
for  obtaining  the  pmf  of  the  random  variable  X  =  Xt  +  X2  is  as  follows. 

Step  1.  Transform  <p(x\,  x2)  to  -k{x\,x)  by  a  substitution  of  variables  (x2  =  x  —  xf 
so  that 


7r(xi,x)  =  (j>(x  i,x  —  xf). 


(29) 
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Step  2.  Sum  n{x\,x)  over  all  Xi  to  obtain  h(x),  the  pmf  for  X  by 

h(x)  =  P{X  =  x}  =  ir{xi,x).  (30) 

Xl 

Let  Sx  =  {x. :  x  =  Xi  +  x2j  where  x  =  (xi,x2)  €  R2.  The  probability  that  X 
assumes  the  value  x  in  its  sample  space  is  given  by 

h(x)  =  P{X=x}=J2^)-  (31) 

xs  sx 

Using  this  approach,  the  distribution  of  the  sum  of  two  random  variables  can  easily 
be  approximated,  even  if  the  component  variables  are  not  independent.  The  pmf 
approximation  is  obtained  from  the  multivariate  kernel  estimate  by  converting  the 
joint  densities  into  point  masses  over  a  large  (but  finite)  set  of  points  as  described  by 
(20)  in  the  univariate  case. 

5.1.  Validation  of  univariate  estimates 

To  demonstrate  the  viability  of  the  univariate  kernel  estimator  in  approximating 
the  component  characteristic  distributions  and  gh  i  =  1,  2,  samples  of  100  observa¬ 
tions  were  taken  from  normal,  symmetric  triangular  and  uniform  populations  (par¬ 
ameters  are  given  in  table  1).  These  populations  were  generated  using  the  inverse 
transform  method.  The  Kolmogorov-Smirnov  (K-S)  goodness-of-fit  test  uses  the 
maximum  absolute  deviation  between  the  estimated  and  true  cumulative  distribution 
functions  (CDFs)  for  a  fixed  number  of  points  (m)  as  the  test  statistic.  If  the  calcu¬ 
lated  test  statistic  exceeds  the  critical  value,  then  the  null  hypothesis  that  the  dis¬ 
tributions  are  equal  is  rejected.  In  particular,  for  m  >  40,  the  K-S  critical  value  is  (cf. 
Conover,  1980) 

9  =  1.36m-1/2.  (32) 

In  the  univariate  experiments,  the  kernel  and  actual  CDFs  were  compared  at  151 
points  which,  at  the  0.05  level,  yielded  a  maximum  allowable  deviation  in  probability 
of  9  =  0.111.  For  each  distribution  type,  five  random  samples  (of  100  observations 
each)  were  drawn  and  kernel  estimates  constructed  for  each  case.  Table  1  sum¬ 
marizes  the  cases  and  reports  the  calculated  K-S  test  statistic.  It  is  seen  that,  in 
each  case,  we  fail  to  reject  the  hypothesis  that  the  distributions  are  equal. 

Figure  2  depicts  the  distribution  function  of  the  kernel  estimate  and  the  actual 
CDF  of  a  symmetric  triangular  random  variable  with  minimum  value  of  50,  mode  of 
75  and  maximum  of  100.  With  a  maximum  absolute  deviation  of  0.041 1  in  prob¬ 
ability,  it  is  seen  that  the  kernel  estimate  closely  matches  the  true  triangular  distri¬ 
bution  function. 


Sample  no. 

Normal  (75,  69.44) 

Triangular  (50,  75,  1000) 

Uniform  (50,  100) 

1 

0.0300 

2 

0.0638 

3 

0.0486 

0.0634 

0.0298 

4 

0.0181 

0.0425 

0.0801 

5 

0.0824 

0.0497 

0.0404 

Table  1.  Maximum  absolute  deviation  in  probability  of  actual  versus  kernel  CDFs. 


Statistical  tolerance  analysis  for  component  characteristics 


347 


Figure  2.  Sample  distribution  function  X  ~  7750,  75,  100). 


It  was  further  desired  to  test  the  ability  of  the  univariate  CDFs  to  predict  the 
appropriate  g  (gj,g2)  as  Per  the  technique  presented  in  section  3.  However,  in  each 
case,  the  true  population  standard  deviation  was  replaced  by  the  sample  standard 
deviation,  which  is  based  on  100  observations.  In  this  comparison,  the  approximated 
g  was  compared  with  the  known  theoretical  g  for  the  normal,  triangular  and  uniform 
distributions  assuming  C*  =  1.  It  is  known  that  for  the  normal  distribution,  g  =  6.0 
will  cover  99.73%  of  the  distribution  with  0.135%  on  each  tail.  For  the  triangular 
distribution,  g  =  4.899  will  yield  100%  coverage  of  the  distribution,  while  it  is  3.464 
for  the  uniform  distribution  with  100%  coverage.  Table  2  summarizes  the  results 
from  the  univariate  experiment. 

5.2.  Non-normal  component  characteristics 

To  demonstrate  the  efficacy  of  the  kernel  estimator  for  determining  the  distri¬ 
bution  of  the  assembly  characteristic  X  when  the  component  characteristics  are  non¬ 
normal  and  independent,  three  distributions  were  assumed  for  the  component  char¬ 
acteristics  and  all  combinations  tested.  A  summary  of  the  distribution  parameters  is 
given  in  table  3. 

The  cumulative  distribution  function  of  the  assembly  characteristic  was  obtained 
via  Monte  Carlo  simulation  (population  size  12  000).  The  K-S  one-sample  test  was 
again  applied  to  assess  statistical  equality  of  the  simulated  and  kernel  distributions. 
Cumulative  distribution  functions  were  compared  at  50  distinct  points  which,  by 


Sample  no. 

Normal  (99.73%) 

Triangular  (100%) 

Uniform  (100%) 

Estimate 

Actual 

Estimate 

Actual 

Estimate 

Actual 

1 

5.777 

6.000 

5.172 

4.899 

3.815 

3.464 

2 

5.933 

6.000 

4.804 

4.899 

3.779 

3.464 

3 

6.231 

6.000 

5.047 

4.899 

4.986 

3.464 

4 

6.071 

6.000 

4.713 

4.899 

3.870 

3.464 

5 

5.644 

6.000 

4.738 

4.899 

3.749 

3.464 

Table  2.  Estimated  versus  actual  g  for  component  characteristics. 


348 


J.  P.  Kharoufeh  and  M.  J.  Chandra 


Case  no. 

X\  distribution 

X2  distribution 

1 

triangular  (50,  75,  100) 

uniform  (50,  100) 

2 

triangular  (50,  75,  100) 

exponential  (75_  ) 

3 

uniform  (50,  100) 

exponential  (75_1) 

Table  3.  Component  characteristic  distribution  parameters. 


Assembly  characteristic 


Figure  3.  Sample  distribution  function  estimate  with  100  observations. 


(32),  gives  a  critical  value  of  9  =  0.192  for  the  K-S  test.  Figure  3  shows  a  typical 
example  of  the  performance  of  the  kernel  technique  in  estimating  the  distribution 
function  of  the  assembly  characteristic  when  the  first  component  is  symmetric  trian¬ 
gular  and  the  second  is  exponentially  distributed. 

Table  4  summarizes  the  results  of  all  goodness-of-fit  tests  for  10  trials  under  each 
case  and  it  is  observed  that,  in  every  instance,  we  fail  to  reject  the  null  hypothesis  that 
the  estimated  distribution  is  equal  to  the  true  distribution  at  the  0.05  level.  Hence,  we 
may  be  confident  that  the  distribution  functions  used  for  the  tolerance  allocation 
technique  will  yield  reliable  estimates.  The  next  section  considers  the  case  where  the 


Sample  no. 

Triangular  +  uniform 

Triangular  +  exponential 

Uniform  +  exponential 

1 

0.1264 

0.0818 

2 

0.0789 

0.0796 

3 

0.0692 

0.0611 

0.1074 

4 

0.0791 

0.1068 

5 

0.1203 

0.0826 

6 

0.1264 

0.1364 

7 

0.1387 

0.1082 

8 

0.0847 

0.1017 

0.1175 

9 

0.0775 

0.0998 

0.1137 

10 

0.0798 

0.0906 

0.1349 

Table  4.  Maximum  absolute  deviation  in  probability,  simulated  versus  estimated  CDF.  The 
critical  value  is  0.192  at  the  0.05  level. 
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two  component  characteristics  are  normally  distributed  and  correlated  random 
variables. 

5.3.  Correlated,  normal  component  characteristics 

In  this  experiment,  we  attempt  to  characterize  the  CDF  of  the  assembly  char¬ 
acteristic,  X  when  the  two  component  characteristics  are  identically  distributed 
normal  random  variables  which  are  correlated.  Recall  that  the  sum  of  two  correlated 
random  variables  is  again  normal  with  mean  equal  to  the  sum  of  the  individual 
means  and  variance  given  by  (6).  We  again  compare  our  discretized  kernel  estimates 
with  the  true  distribution  of  X1  +  X2  by  means  of  the  K-S  test.  In  each  case,  the 
kernel  estimate  of  the  joint  density  function  was  constructed  using  a  random  sample 
of  50  observations  of  (Xl,  X2).  The  normal  random  variables  have  a  mean  of  75  and 
standard  deviation  of  10.  The  resulting  CDFs  were  compared  at  50  distinct  points 
(not  the  observations)  to  test  the  hypothesis  that  the  distributions  are  equal  at  the 
0.05  level.  Table  5  gives  a  summary  of  the  empirical  results. 

Figure  4  depicts  the  estimated  distribution  function  (‘Kernel’)  of  the  assembly 
characteristic  X  versus  the  actual  population  CDF  (‘Actual’).  It  is  important  to  note 
that  each  of  these  was  constructed  using  only  50  sampled  components  from  the 
populations  of  Xl  and  X2  This  is  significant  because  a  relatively  small  sample  size 
is  used  to  characterize  the  distribution,  even  when  the  components  are  correlated. 

6.  Conclusions  and  future  directions 

This  paper  has  presented  an  approximation  technique  for  statistical  tolerancing 
analysis  that  allocates  component  tolerances  for  a  given  desired  assembly  tolerance. 
The  main  advantage  of  the  simple  approach  is  that  it  is  a  data-driven,  non-para- 
metric  approach  that  allows  the  process  data  to  dictate  the  appropriateness  of  com¬ 
ponent  characteristic  distributions  and  the  assembly  characteristic  distribution. 
Thus,  the  typical  assumptions  of  independence  and  normality  can  possibly  be 
relaxed  when  the  number  of  components  being  combined  is  small  (i.e.  k  <  10). 

To  validate  the  technique’s  ability  to  characterize  the  assembly  characteristic 
distribution,  we  examined  two  cases  for  two  components:  (1)  the  two  component 
characteristics  are  independent  and  not  normally  distributed,  and  (2)  the  component 
characteristics  are  normally  distributed  but  highly  correlated.  In  the  first  case,  the 
results  indicate  that  kennel  density  estimates  can  characterize  the  component  distri¬ 
butions  and  the  overall  assembly  characteristic  distribution  with  reasonable  preci¬ 
sion.  The  maximum  absolute  deviations  (in  probability)  are  not  statistically 
significant  at  the  0.05  level  for  three  different  distribution  types.  In  the  second 
case,  it  was  found  that  the  kernel  estimate  performs  well  when  estimating  the  assem¬ 
bly  characteristic  distribution  even  when  the  component  characteristics  are  highly 
correlated  normal  random  variables.  A  logical  extension  is  to  consider  the  case  of  a 
small  number  of  highly  correlated,  non-normal  component  characteristics.  This 
work  is  currently  underway. 

Both  experiments  indicate  that  by  using  actual  process  data,  it  may  be  possible  to 
relax  the  typical  assumptions  made  in  statistical  tolerancing  analysis,  namely  the 
independence  and  normality  assumptions  for  the  component  characteristics.  Of 
course,  as  with  any  statistical  estimator,  there  are  pros  and  cons  for  using  the 
kernel  density  estimates.  One  disadvantage  is  the  process  of  selecting  the  smoothing 
parameter,  /?.  Although  this  process  may  be  automated  to  some  degree,  it  is  instruc¬ 
tive  to  ‘fine  tune’  the  parameter  to  achieve  a  better  fit  of  the  data,  as  a  different 
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Test  no. 

Correlation 

Maximum  deviation  in  probability 

Critical  value 

1 

0.5000 

0.0731 

0.1920 

2 

0.7070 

0.0713 

0.1920 

3 

0.9480 

0.0790 

0.1920 

4 

0.9950 

0.0837 

0.1920 

Table  5.  Test  results  for  X  =  X\  +  X2,  where  X\  and  X2  normal  and  correlated. 


Assembly  Characteristic  Assembly  Characteristic 

Figure  4.  Assembly  characteristic  distribution  when  components  have  a  correlation  coeffi¬ 
cient:  (a)  0.500,  (b)  0.707,  (c)  0.950,  (d)  0.995. 


smoothing  parameter  yields  possibly  dramatically  different  distributions.  One 
obvious  advantage  of  using  the  kernel  estimator  is  that  we  make  no  distributional 
assumptions  with  regard  to  the  component  or  assembly  characteristics,  and  do  not 
lose  much  in  the  discretization  of  the  distributions  so  long  as  the  smoothing  par¬ 
ameter  is  chosen  appropriately.  Hence,  with  some  good  historical  process  data,  it 
may  be  possible  to  determine  quickly  the  appropriate  component  tolerances  to 
achieve  the  overall  desired  assembly  tolerance.  Even  with  a  small  number  of  samples, 
it  is  possible  to  find  reasonable  values  for  the  distribution  functions. 
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As  indicated  here,  the  concept  of  process  capability  ratios  (capability  indices)  is 
very  closely  related  to  the  statistical  allocation  of  component  tolerances.  In  particu¬ 
lar,  it  was  seen  that  the  target  Cp  index  was  used  explicitly  to  find  equations  needed 
to  allocate  the  known  assembly  tolerance  among  k  components.  This  procedure  can 
be  easily  extended  to  incorporate  other  target  capability  indices  such  as  Cpk  or  Cpm 
by  using  the  maximum  deviation  allowable  (,?,-)  between  the  process  mean  and  nom¬ 
inal  dimension  d These  extensions  are  applicable  even  for  non-normal  or  correlated 
component  characteristics  since  the  approach  is  driven  by  the  actual  process  data  to 
compute  the  needed  values,  c  and  ch  as  indicated  in  section  3.2. 

In  the  future,  it  would  be  instructive  to  investigate  more  sophisticated  kernel 
estimators  that  are  adaptive.  In  this  paper,  the  smoothing  parameter  was  fixed  in 
both  coordinate  directions,  but  it  could  be  varied  by  choosing  a  vector  of  smoothing 
parameters  in  each  direction.  Such  an  adaptive  method  may  lead  to  better  estimates 
of  the  joint  density,  and,  hence,  of  the  assembly  characteristic  distribution.  The 
kernel  does  lend  itself  to  extension  to  higher  dimensions,  and  we  are  currently 
exploring  this  to  include  more  than  two  characteristics.  Ultimately,  the  technique 
will  be  incorporated  into  a  cost-minimization  scheme  and  compared  against  existing 
techniques. 

Acknowledgements 

The  authors  acknowledge  the  helpful  comments  of  two  anonymous  referees 
whose  suggestions  helped  to  clarify  several  concepts  in  the  paper. 


References 

Barker,  A.  M.,  1989,  Statistical  tolerancing:  methods  to  improve  quality.  Proceedings  of  the 
43rd  Annual  Quality  Congress  Transactions,  8-10  May,  Vol.  43,  Toronto,  pp.  490M95. 
Bender,  A.,  Jr,  1968,  Statistical  tolerancing  as  it  relates  to  quality  control  and  the  designer. 
SAE  Transactions ,  77,  1965-1971. 

Bjorke,  O.,  1989,  Computer-aided  Tolerancing ,  2nd  edn  (New  York:  ASME  Press). 

Burr,  I.  W.,  1958,  Some  theoretical  and  practical  aspects  of  tolerances  for  mating  parts. 
Industrial  Quality  Control,  15,  18-22. 

Conover,  W.  J.,  1980,  Practical  Nonparametric  Statistics  (New  York:  Wiley). 

D’Errico,  J.  R.  and  Ziano,  N.  A.,  1988,  Statistical  tolerancing  using  a  modification  of 
Taguchi’s  method.  Technometrics,  30,  397-405. 

Evans,  D.  H.,  1974,  Statistical  tolerancing:  the  state  of  the  art.  Part  I.  Journal  of  Quality 
Technology ,  6,  188-195. 

Freund,  J.  and  Walpole,  R.,  1987,  Mathematical  Statistics,  4th  edn  (Englewood  Cliffs: 
Prentice-Hall). 

Fukunaga,  K..,  1972,  Introduction  to  Statistical  Pattern  Recognition  (New  York:  Academic 
Press). 

Mansoor,  E.  M.,  1963,  The  application  of  probability  to  tolerances  used  in  engineering 
design.  Proceedings  of  the  Institution  of  Mechanical  Engineers,  178  ,  29-44. 

Nigam,  S.  D.  and  Turner,  J.  U.,  1995,  Review  of  statistical  approaches  to  tolerance  analysis. 
Computer-Aided  Design ,  27,  6-15. 

O’Connor,  M.  A.  and  Srinivasan,  V.,  1997,  Composing  distribution  function  zones  for 
statistical  tolerance  analysis.  Proceedings  of  the  5th  CIRP  International  Seminar  on 
Computer-Aided  Tolerancing  (CAT),  pp.  13-24. 

Parkinson,  D.  B.,  1982,  The  application  of  reliability  methods  to  tolerancing.  Transactions  of 
the  ASME  Journal  of  Mechanical  Design,  104,  612-618. 

Parkinson,  D.  B.,  1984,  Tolerancing  of  component  dimension  in  CAD.  Computer-Aided 
Design ,  16,  25-31. 

Scott,  D.  W.,  1992,  Multivariate  Density  Estimation,  Theory,  Practice,  and  Visualization 
(New  York:  Wiley). 


352 


Statistical  tolerance  analysis  for  component  characteristics 


Silverman,  B.  W.,  1986,  Density  Estimation  for  Statistics  and  Data  Analysis  (New  York: 
Chapman  &  Hall). 

Zhang,  H.  C.  and  Huq,  M.  E.,  1992,  Tolerancing  techniques:  the  state-of-the-art 
International  Journal  of  Production  Research ,  30,  2111-2135. 

Zhang,  C.,  Luo,  J.  and  Wang,  B.,  1999,  Statistical  tolerance  synthesis  using  distribution 
function  zones.  International  Journal  of  Production  Research,  37,  3995^4006. 


